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ABSTRACT 



Electromagnetic scattering from trees and vegetation is of prime importance in 
radar and remote sensing. The actual problem of scattering from trees is rather 
complicated and involves three dimensional scattering from lossy, electrically large, 
and randomly oriented objects. 

In this thesis, the radar cross section of a planar fractal tree is considered. 
Although a planar tree is far from being real, scattering from it sheds light on the 
scattering phenomenon from an actual tree. The planar tree is generated using 
fractal geometry and its branches are considered perfectly conducting. The tree is 
illuminated by a plane wave and the problem is solved using the moment method. 
Data is presented for the radar cross section for different branching angles of the 
tree and at different frequencies. 
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I. INTRODUCTION 



A. NEED FOR THE STUDY 

The trees existing in the natural world are fractal anisotropic. They are made 
up of long, intersecting and lossy objects. The geometry of these objects is not easy 
to set up as they are randomly oriented in a three dimensional space. The use of 
fractals facilitates the modeling of semi— randomly distributed structures. 
Mandelbrodt [Ref. 1] has shown that a number of naturally occurring phenomenon 
such as coastlines, clouds, trees, etc. are fractal in nature. For instance, when a 
branch is divided into two (or more), the ratio of the length of the subbranches to 
the main branch length remains constant. Furthermore, the branching angle also 
remains the same. 

In this thesis, scattering from planar fractal trees is considered. The radar 
cross section of a real fractal tree is a complicated scattering problem. The analysis 
of this problem in a two dimensional space gives an approximate idea of the radar 
cross section of a real tree. 

The radar cross section of an object is a quantitative measure of the ratio of 
the power density that is received and scattered by the object to the power density 
of the electromagnetic wave that illuminates that object. The radar cross section is 
independent of the range of the object for the far— field situation. The theoretical 
definition of the radar cross section 'V' is given by the formula: 

a = 47tR^ lim 1—^—1 
R - oo 1 fi 1 1 



1 



where E 1 is the incident electric field vector, fi s is the scattered electric field vector, 
and R is the distance between the scattered object and the point of observation. The 
radar cross section has dimensions of area. Usually, it is expressed in square 
wavelengths. 

B. STATEMENT OF THE PROBLEM 

This thesis investigates the radar cross section of planar fractal trees. These 
trees are composed of planar thin strip dipoles of arbitrarily dimensions and 
orientations in a plane. These structures are excited by plane waves of various 
frequencies. 

The first step in solving the problem is to calculate the induced current 
distribution on each strip. The calculation of the current distribution is based on the 
theory of the moment method and requires a knowledge of the impedance between 
any two of these strips as well as the voltage on each planar strip due to the 
incident electric field. The basic concepts and the calculation of the current 
distribution are described in Chapter 2. 

In Chapter 3, the development of a FORTRAN program, is presented. The 
evaluation of the radar cross section requires the knowledge of the scattered electric 
field due to the induced currents on the planar strips. The program computes the 
scattered electric field and then the radar cross section of that structure. The details 
of these calculations are also presented in this Chapter. 

The computer models that are used by the developed program are presented in 
Chapter 4. Their generation is based on the fractal geometry. An existing and 
modified program is used to generate the geometry of the planar fractal trees in 
order to be used as input in the developed program. 



2 



The numerical results of the radar cross section of a single planar dipole and a 
a number of planar fractal trees are presented in Chapter 5. The scattering from a 
single dipole is compared with standard results for a similar case. The limitations of 
the developed program are also presented in this Chapter. 

In Chapter 5, the conclusions of the radar cross section results and 
recommendations are presented. The two programs that are used for this 
investigation are listed in the Appendices. 
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II. METHOD OF MOMENTS THEORY 



In this section of the study, the basic concepts of the moment method theory 
are presented. This theory is used in the development of the RCS program to find 
the current distribution on a planar strip due to an incident plane wave. 

For a given structure consisting of planar dipoles the impedance between any 
two of them is calculated from the knowledge of the geometry and the wavelength of 
the incident plane wave. The voltage on each dipole is calculated from the 
knowledge of the characteristics of the incident electric field. The induced current 
distribution on each dipole of the structure is determined from the calculated 
impedance and voltage using the method of moments theory. 

A. GENERAL THEORY 

The method of moments is a numerical procedure for solving integral 
equations of the form: 

f b 

f(x')K(x,x')dx' = g(x), a<x<b (eqn2.1) 

Ja 

where f(x') is an unknown function, K(x,x') is a known Kernel or Green's function, 
and g(x) is a given function. This procedure reduces the integral equation (eqn 2.1) 
to a system of simultaneous linear algebraic equations in terms of some unknown 
coefficients. This method requires that the function f(x') be approximated by a 
series of N expansion functions or " basis functions ", such that 



4 



n=l,2,.. 



,N 



(eqn 2.2) 



f(x')x Ea n f n (x'), 

n = l 



where the domain of f n (x') is the same as that of f(x') and a n 's are the complex 
unknown expansion coefficients. 

There are two types of basis functions. The subdomain functions, which are 
nonzero over a part of the domain of the unknown function f(x'), and entire domain 
functions being nonzero over the entire domain of f(x'). In antennas, some 
commonly employed subdomain basis functions are the piecewise sinusoid functions, 
the unit height— pulse functions and the piecewise triangular functions. 

The subdomain procedure requires subdivision of the structure into N 
nonoverlapping segments. Figure 2.1 shows a segmented line where the segments are 
assumed to be collinear and of equal length, although this condition is not necessary. 



S-e INTERVAL 



a Xi X2 x 3 X4 x 5 b 



Figure 2.1 Segmented Line. [From Ref. 2] 



Figure 2.2 shows a subdomain unit height— pulse function which produces a 
staircase representation of the unknown function f(x'). Figure 2.3 shows a 
subdomain sinusoid basis function and the representation of the function f(x'). 
Figure 2.4 shows a subdomain triangular basis function producing a smoother 
representation of the function f(x') than the case of the unit height-pulse basis 
function. 

The use of entire— domain basis functions does not require any segmentation of 
the structure. One of the most most commonly used basis functions of this kind is 
the sinusoidal basis functions. 
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Figure 2.2 Unit Height-Pulse Basis Function and a Staircase 
Representation of f(x'). [From Ref. 2] 





Figure 2.3 Sinusoid Basis Function and a Representation 
off(x / ). [From Ref. 3] 





Figure 2.4 Triangular Basis Function and a Representation 
of f(x'). [From Ref. 2] 
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The substitution of the function f(x') by a sequence of N basis functions leads 
to one equation of N unknowns of the expansion coefficients a n , which can be found 
by using N linearly independent equations. These equations are set up by the use of 
"testing or weighting" functions u„,(x), m=l,2,....N. At this point the definition of 
the inner product is required. The inner product (f, g) between any two functions or 
vectors f, g is a scalar operation defined as 



where S is the surface of the structure that is analyzed [Ref 4]. The inner product of 
the selected testing functions w m (x), with the two sides of the original integral 
equation leads to the equation: 




(eqn 2.3) 




(eqn 2.4) 



or 




& n = l 



m = 1,2, ,N 



(eqn 2.5) 



The use of the above definition for the inner product, yields: 



The following substitutions 



f b 

Y m = g(x)u; m (x)dx, m=l,2 ,N 



(eqn 2.7) 



a 



and 



r b r b 

Zmn = a^(x)dx f n (x' )K(x,x ' )dx' , m,n=l,2,....N 



a J a 



(eqn 2.8) 



result in a matrix equation of the form: 



[Z mn ] * [&n] — [^m] 



(eqn 2.9) 



The unknowns [a n ] can be obtained through a matrix inversion: 




(eqn 2.10) 
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where [Z mn ] is the inverse matrix. The column vector [V m ] depends upon the given 
function g(x) and the selected testing functions w m (x). The matrix [Z mn ] depends 
upon the known kernel K(x,x') and both the selected basis and testing functions. 
Once the expansion coefficients are known, the function f(x') is also known. 

The choice of basis and testing or weighting functions is based upon experience 
and the rule is that their number has to be the same. The procedure of using the 
same basis and testing functions is called the " Garlekin's method ". 

B. APPLICATIONS TO EM THEORY 

A number of problems in electromagnetic radiation and scattering can be 
solved by the method of moments. In the present case, the simple case of a perfectly 
conducting object 11 situated in free space " is considered. 

The basic problem is to investigate the case when an object is illuminated by 
fields of known impressed electric and magnetic currents (J 1 , M 1 ). In the absence of 
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the object, the impressed currents radiate the assumed known incident electric and 
magnetic fields (E 1 , H 1 ). In the presence of this object, the impressed currents 
radiate the unknown total fields (e\ H^). 

The integral equation is obtained by the surface equivalence principle, 
replacing the object by free space together with the electric surface current density 
J = nxH t (eqn 2.11) 

where J exists on the entire surface S of the object and n is the unit vector normal 
to that surface. In free space, J radiates the scattered fields E s , H s 

E s = E t — E 1 (eqn 2.12) 

II s = H t — E 1 (eqn 2.13) 

The boundary conditions for this case enforce the total tangential electric field 
on the surface S to zero: 

n x ( E s + E 1 ) = 0 (eqn 2.14) 

This is an integral equation for J since the scattered electric field E s can be written 
as an integral over S of the dot product of J and the dyadic free space Green's 
function. [Ref. 5] 

As the geometry of the object is known, it is more convenient to use the total 
current I instead of the total current density J. So, the problem is to find the 
unknown current I induced by the incident electric field. The method of moments 
solve these kind of problems by the following procedure. 

The first step is to expand the unknown current I in terms of some basis set: 

N 

Is EI n F n (eqn 2.15) 

n = l 

where the I n are the sequence of N unknown complex coefficients, and the F n is a 
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sequence of N known modes or basis functions. The best choice of F n for a given 
problem could be quite involved and is discussed in [Ref. 4, pp 308—310]. 

The second step is to select the testing or weighting functions u; m , m=l,2,....N. 
These can be identical with the basis functions or different. The inner product of the 
sequence of N weighting functions u^, with both sides of the integral equation gives 
a N x N system of simultaneous linear algebraic equations of the symbolic form: 



where I is the current column vector whose N components give the values of I n 
and [Z] is the N x N impedance matrix given by the equation 



where S is the surface of the structure being analyzed. The impedance matrix [Z] is 
always symmetric, and, for the special case of a thin dipole instead of an arbitrary 
object, is also a Toeplitz matrix. In a Toeplitz matrix, Z mn depends only on | m— n | . 

Generally, [Z] is dependent only on the geometry and material composition of 
the scatterer, but not on the incident fields [Ref 5]. The right-hand side of the last 
equation is the voltage vector whose N components give the corresponding mode 
voltage. The voltage vector depends only on the excitation, i.e., the incident electric 
field. The dimensions of the elements of [Z] and [V] are volt-amps (VA), while the 
elements of [I] are dimensionless. 

The solution of the last matrix equation is the column vector [I], whose 
elements represent the complex coefficients I n . As the total current I was expanded 
by a sequence of N known modes or basis functions and the complex coefficients I n 
are known, the total current and so the total current density J is also known. 



[Z].(I] = (V] 



(eqn 2.16) 




(eqn 2.17) 
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Although the choice of weighting functions is free, it has to be considered that 

2 

the matrix equation being solved requires the evaluation of N terms and each term 
requires two or more integrations. When these integrations are to be done 
numerically, the computations become complicated. There is a way to reduce this 
complexity by choosing as weighting functions the Dirac delta functions. This is the 
method of point-matching in which delta functions are enforced only at discrete 
points on the surface S. The results of this method can be quite accurate especially 
when the discrete points are selected to be equally spaced. The solution satisfies the 
electromagnetic boundary conditions (e.g., vanishing tangential electric fields on the 
surface of an electric conductor) only at discrete points [Ref 4]. Between these 
points the boundary conditions may not be satisfied. In this case it is required to 
define the deviation as a residual (e.g., residual=AE| ta n = E s | tan + E'ltan 4 0 on 
the surface S of an electric conductor) and use the method of weighted residuals so 
that the boundary conditions will be satisfied in an average sense over the entire 
surface S. 
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III. ANALYSIS AND DEVELOPMENT OF RCS PROGRAM 



In this chapter, the basic steps that are involved in the calculation of the 
Radar Cross Section (RCS) of a planar structure composed of conducting strips are 
presented. As mentioned earlier, the fractal tree is modeled as consisting of planar 
thin strips. The radar cross section of the fractal tree is calculated using the method 
of moments. A Fortran program is developed to calculate the RCS of the tree for a 
specified geometry and at a given frequency. 

The moment method discretizes an integral equation to a matrix equation of 
the form: 

[Z]-[I] = [V] (eqn 3.1) 

where [Z] is the impedance matrix whose elements represent the mutual impedance 
between any two dipoles of a given structure depending upon the wavelength and 
geometry of the structure, [V] is the voltage matrix w T hose elements correspond to 
the voltage on each dipole due to excitation ( incident electric field ), and [I] is the 
unknown matrix whose elements represent the induced current on each dipole. 

A. DEVELOPMENT OF RCS PROGRAM 

The numerical results for RCS are obtained from a Fortran program that 
calculates backscattered RCS of a planar structure consisting of arbitrarily oriented 
dipoles. For convenience, the structure will be assumed to be in the y-z plane of an 
xyz cartesian coordinate system. 

Figure 3.1 shows the structure to be investigated. A large number of planar 
dipoles of variable lengths, widths, and orientations in the y— z plane is 



15 



illuminated by a plane wave with electric field linearly polarized and characterized 
by the angles ^>o and 9q. It is required to calculate the backscattered RCS from this 
structure. 

The incident electric field E 1 is given by the formula 



P = e. e - j( kx ‘ x + Vy + kz ' z ) 


(eqn 3.2) 


where £ is the propagation vector and 




k x ^+ ky^+ k z ^ = ko^ = u?hq€o 


(eqn 3.3) 


CO 1 

II 

o 


(eqn 3.4) 


k x = kosin#ocos</?o 


(eqn 3.5) 


k y = kosin^osin^o 


(eqn 3.6) 


k z = k o cos0 o 


(eqn 3.7) 


e = E x x + E y y + E z z 


(eqn 3.8) 


The electric field E 1 lies on the plane perpendicular to 


the direction of 


propagation of the plane wave. Hence, e-k = 0 implies that 




E x k x + Eyky 4- E z k z = 0 


(eqn 3.9) 


The substitution of E y and E z by the variables a and b (independent variables) gives 


the coordinate E x as a dependent variable 




E x _ (ak y + bk z ) 
k x 


(eqn 3.10) 


Equation 3.1 the becomes 




ff-Ly + bz (...aky + bk z ) - ] e ~j( k,x + k y y + k 2 z ) 

1 k x J 


(eqn 3.11) 
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Figure 3.1 Geometry of RCS Program. 
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The RCS program that has been developed makes the following assumptions: 

1. The dipoles are very thin planar strips. The width of the dipoles is 
assumed to be electrically small so that only the axial current is 
significant. Further, there is no current variation along the width of 
the dipoles. 

2. The dipoles are not intersecting. 

3. The dipoles are perfect conductors. 

The program uses the theory of moment method and the way that it solves the 
problem can be understood if a single planar thin dipole in the y— z plane is 
considered. Each dipole is subdivided into equal segments. Overlapping Piecewise 
Sinusoidal (PWS) modes are assumed to exist on the dipole. The length of each 
PWS mode is equal to the length of two segments. Figure 3.2 shows the m t h PWS 
mode. The PWS mode has a length 2h m , and a width 2w„). The coordinates of its 
center are (y^, z m ), and it is oriented at an angle ij> m measured from the z axis. 

The incident electric field induces current along the axis ( of that mode. In 
Figure 3.2, the induced current J varies sinusoidally along the axis of the m t h mode. 
A new coordinate system t/—( is introduced, where ( is the axis of the dipole and rj is 
the axis perpendicular to (. The coordinate system is obtained by rotating the 
y-z system about the x— axis through an angle 90°— Figure 3.3 shows the details 
of this rotation. 

The relation between the coordinates of the center of the mode along the axis 
of the two orthogonal systems is found to be: 

y = y m + Csin(^m) - ^cos(^ m ) (eqn 3.12) 

z = z m + Ccos(^m) + rj sin(^,n) (eqn 3.13) 
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The corresponding vector equation is: 
y = Csin(^m) - rj cos(^m) 
z = Ccos(^m) + ^sin(^m) 



(eqn 3.14) 
(eqn 3.15) 




Figure 3.2 Geometry of the m t h PWS Mode of the Structure. 
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Figure 3.3 Transformation of PWS Mode's Center Coordinates. 
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L Voltage Equations 

In this section, expressions for the elements of the voltage matrix [V] 



due to the incident electric field, are presented. Each element of this matrix 
corresponds to a PWS mode of the structure that is investigated. The size of the 
voltage matrix is equal to the total number of modes of the structure. 

The induced current density on the m t h mode of the structure, due to 
incident electric field, is 

J = -i- _ 1 C l )) (eqn 3.16) 

2w m sin(k 0 h m ) 

where 2w m is the width and 2h m is the length of the m t h mode. For each mode, the 
induced current variation will be sinusoidal along its axis, being maximum at the 
center and zero at the end points. The corresponding voltage V m is 



where t m is the angle between the z axis (reference for angle measurements) and the 
axis of the mode and y m , z m are the coordinates of the center of the mode along the 




(eqn 3.17) 



From equations 3.11 and 3.16, it is seen that 



ttj = [a sin(Vb) + b cos(«] e"j[ k ^" + + k *< z " +<“*(*>)] 



(eqn 3.18) 
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y and z axis respectively. The substitution into equation 3.17, gives the following 
form of V m : 



Q— j(kyym + k z Z m ) p h m 

= sin( k0 h„) [aSinW " ) + b C0SW>)1 j-h. sin ( k o(hm— f c] ))dC 

(eqn 3.19) 

where k<- = k y sin(V>m) + k z cos(^m). A closed form evaluation of the integral in this 
equation leads to final form of V m 



v, = f° v °" 2 [ cos(koh„) - cos^h.) ] 


k; + ± ko 


k; - ko 


(eqn 3.20) 


Vm = V om hm sin(kohm) 


k ; = ± k 0 
(eqn 3.21) 


where 




e ~ j(^yym + k z z m ) p 

Vom = a Sin(V>m) + b cos(^m) 

sin(k 0 h m ) L J 


(eqn 3.22) 



2. Radiation Equations 

In this section, expressions for the far-zone scattered fields due to the 
m t h PWS mode are developed. For far— field observations the electric field is given 
in spherical coordinates by the following equations [Ref 4]: 
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E * 0 



(eqn 3.23) 



E# — + , N ) 

4 7TT ^ 

E , + tei jkor - 



4 7TT 



-( L^+ ^N^,) 



where 



(eqn 3.24) 
(eqn 3.25) 



= JJ [M x cos#eos(^ + M y cosfciny>-M z sin0] e + ^° r cos ^ m ds' 

(eqn 3.26) 

= Jj^ [-M x sinv? + Mycosp] e +j k ° r ' cos V’m dg / 

(eqn 3.27) 

= JJ [J x cos#cos<^ + Jycosfein^ - J z sin0] e + ^° r cos ^ m ds' 



(eqn 3.28) 

N^= J| [-J x sint^ + Jycos^] e + ^° r cos ^ m ds' 

(eqn 3.29) 

= 120- (eqn 3.30) 

The quantities J x , J y , J z , are the components of the electric current 
density J s that are induced on the m t h mode over the surface S, and the quantities 
M x , M y , M z , are the coordinates of the magnetic current density M$ over the surface 
S. As the structure is on the y-z plane and M$ is zero, the quantities J x , L^, are 
zero. Equations 3.24 and 3.25 become: 



|8S ~jk 0 e j k ° r i? N 
4 7TT 



4/rr 





(eqn 3.31) 


N 

V 5 


(eqn 3.32) 
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As the current density J m is along the axis, all quantities within the 
integrals for J and M will be expressed in terms of the i system. The surface 
element that is used in all integrals is ds' = dydz = d(dr) = d£, as the width of the 
mode is assumed to be very small in terms of its length. The transformations from 
the y— z system to the (—77 system are made by using the following substitutions: 



r'cos(^ m ) 3 ysinflsin^ 4- zcostf 



where 



and 



y = ym + Csin(^m) 
z = z m + Ccos(^m) 

Jm = Jyy + J Z Z = ( 



2w m sin(k 0 h m ) 



- sin(k 0 (h m — | Cl ) 



where 



J y = sin(^ m ) 

Jz = J £ COS(^ m ) 



(eqn 3.33) 

(eqn 3.34) 
(eqn 3.35) 



(eqn 3.36) 

(eqn 3.37^ 
(eqn 3.38) 



These substitutions and algebraic manipulations give the quantities and 
for the m t h mode in the (- 7 ? system: 



N^m = N^ o 2 2k ° 2 [ C0S ( E mhm) - cos(k 0 h m ) j 

Em — ko 

h m sin(kohm) 



E m / ± ko 

(eqn 3.39) 

Em = ± ko 

(eqn 3.40) 



and 
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Em i- ±ko 




(eqn 3.41) 




E m = ±k 0 



(eqn 3.42) 



where 



^ _ sin^cosflsin^ - cos^mSinfl j ^[koCyniSin^sinvj+ZmCOs^)] 



sin(k 0 hm) 



(eqn 3.43) 




^ _ sin^mcosq? ^[kotymsin&in^+zmcos#)] 



(eqn 3.44) 



E m = ko(sinV , mSin0sinp 4- cosmos#) 

(eqn 3.45) 

In this thesis, only the monostatic radar cross section will be 
considered. In the radiation equations the angles 6 and p represent the orientation 
angles of the scattered electric field E s . Their relation with the incident angles Oq 
and po for the monostatic case is: 

0=tt—9q (eqn 3.46) 

p = 7r — tpo (eqn 3.47) 



If M denotes the total number of PWS modes in the structure under 



investigation, then 



M 




E N 



(eqn 3.48) 



m=l 



M 




(eqn 3.49) 
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The final expression for the equations 3.24 and 3.25 is 



M 



E «= c m E , N <ta 

m=l 

M 

E = C £ N 

* m =l ¥** 


(eqn 3.50) 
(eqn 3.51) 


where 




c _-ik 0 e _ ^ k ° r 7 1 
4 TIT 


(eqn 3.52) 


3. RCS Equations 




When an incident plane wave with electric field E 1 


strikes the object 


and E s is the scattered electric field, the radar cross section a is defined as 




(eqn 3.53) 


r -> °° . 

|E'| 




For the scattered electric field E s , its spherical coordinates E^, and E^ have been 
calculated by the equations 3.50 and 3.51. The incident electric field E 1 is known by 


means of equation 3.11. 




| fi* | 2 - |a| 2 + | b | 2 + |.- ak r.+ bk * | 2 
kx 


(eqn 3.54) 


2 2 2 2 2 2 
|E S | = |E„| + |E^| = | C | [[|N e | + |N^| ] 


(eqn 3.55) 
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The combination of equations 3.11, 3.50-3.52, 3.54-3.55, and 3.53 leads to the final 
formula for RCS: 

a = ^\ IN/+ | N ,„| 2 1 tt , i 

aky + bk z 

k x 

(eqn 3.56) 

This is the formula that the program uses to compute the RCS of a 
planar structure for a given set of incident angles 0o, and <^o- 

B. INPUT-OUTPUT OF RCS PROGRAM 

The RCS program, which is described in Appendix A, is a FORTRAN 
program having two input files. The first input file, INPUT 1, contains the data that 
characterize the incident plane wave and the data that describe the geometry of the 
planar structure whose broadside RCS is measured. The second input file, INPUT, 
contains the set of incident angles Oq, <A)- 

The program reads from the INPUT1 file the following input data: 

1. frequency of the incident plane wave in GHz, 

2. parameters a and b that characterize the polarization of the incident 
electric field. 

3. number of dipoles that the structure consists of, 

4. half length of each dipole in cm, 

5. half width of each dipole in cm, 

6. coordinates in cm of the center of each dipole along the two axes of 
the orthogonal system that the structure lies, 



4x 



|a| + |b| 
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7. orientation angle of each dipole, measured from the vertical axis, 
positive in the clockwise direction, and 

8. number of segments of each dipole. 

As the program reads these input data it generates the geometry of the PWS 
modes, calculates the impedance elements between any two modes of the structure, 
and fills the impedance matrix [Z]. The size of the impedance matrix depends on the 
total number of modes. When this matrix is calculated and filled, the program 
inverts it to [Z]~ * and stores it for later use. 

The next step is to read from the INPUT file the sets of incident angles 9 o, <po 
and for each one it calculates the voltage on each mode, filling the voltage matrix 
[V]. This matrix is a column vector which depends upon the excitation only. Then it 
multiplies the stored inverted impedance matrix [Z] — * by the voltage matrix. This 
yields the current vector [I]: 

[I] = [Z— '-[V] (eqn 3.57) 

As the induced currents are known, the program uses the previously described 

radiation equations and calculates the RCS of the structure corresponding to the 

given set of incident angles 9o, <po ■ The output RCS is normalized to the square of 
o 

the wavelength A , and is given in dB. The program reads the next set of incident 
angles 9q , <po and repeats the same procedure to compute the RCS of the new set. 

Although the input lengths and widths are in cm, the program considers them 
normalized to the wavelength. For accurate results at least 4 segments per 
wavelength are chosen. The selection of the width of each dipole is arbitrarily taken 
as L/W = 33. 
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IV. FRACTAL TREE GENERATION 



The problem of measuring the radar cross section of a real natural tree is 
complicated as it requires an investigation in three dimensional space with very 
long, intersecting elongated objects that are randomly oriented. 

For simplicity the radar cross section (RCS) problem is solved in a two 
dimensional space. RCS is calculated from a planar fractal tree whose branches are 
considered as thin planar dipoles of different lengths and widths. In an actual tree, 
the branches are lossy and in general anisotropic. However, in the present model, 
the tree is comprised of perfectly conducting branches. The geometry of this tree is 
based upon the fractal geometry. 

A. DEFINITIONS 

Fractal is a mathematical set or object whose form is extremely irregular or 
fragmented at all scales [Ref. 6]. The requirement to describe the shape of many 
objects that appear in the natural world, such as trees, mountains, coastlines, etc., 
led to the generation of fractal geometry. As Euclidean geometry cannot give 
mathematical expressions to describe fractal objects, fractal geometry is used to 
describe mathematically many natural patterns. 
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The basic characteristics of fractal objects are [Ref. 7]: 

1. A large degree of heterogeneity. 

2. A self-similar structure over many size scales. Self-similarity refers 
to the general preservation of form or characteristic regardless of the 
scale of observation. 

3. The lack of a well-defined (characteristic) scale. 

The geometric characteristics of fractal objects are useful for describing phenomena 
of nature, such as scattering of objects like landscapes or surface cracks. Although 
these objects are irregular and randomly oriented in nature, they show structural 
similarities on several different discrete size scales [Ref. 7]. 

One measure of structural complexity is the fractal dimension Df. There are 
several definitions of Df depending upon the particular application. For the case of 
fractal trees the fractal dimension Df is the measure of the space that a 

self-similar structure flits, and it varies with the branching levels that the structure 
consists of. 

The most useful terms from fractal geometry that are used to describe a planar 
fractal tree are the following: 

1. The number of branch segments N formed from each preceding branch 
segment. 

2. The constant similarity ratio "r" that relates the fractional reduction 
in segment length for each segment to previous level, this factor is less 
than 1.0. 

In this case, the fractal dimension Df is given by the formula [Ref. 7]: 

^ p r — logftotal branch length) _ log(rN) 
log(aver age branch length) log(l/r) 
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This equation is accurate for symmetric structures or asymmetric structures having 
a normal distribution function [Ref. 7]. More details about fractal geometry are 
discussed in [Ref. 1] and [Ref. 7]. 

B. FRACTAL TREE GENERATION 

In the planar fractal structures that are used for calculating the radar cross 
section, the number of branch segments N formed from each preceding branch 
segment is 2. The principal properties of these trees are that the branches are not 
overlapping, and the angle 9 between any two branches is the same. The 
nonoverlapping between the branches of each fractal model is achieved by choosing 
the branch angle 9 for a given reduction factor by the following empirical formula 
[Ref. 7] 



Minimum 9 = 32.34 * r^'^ (Radians) 

This angle is given in radians and r is the desired reduction factor which is a 
constant for the structure. As the reduction factor r increases the tree fills more 
space. 

Five types of planar fractal trees are considered in this thesis. Each type 
corresponds to a set of values for reduction factor r and minimum branch angle 0 . 
Table 4.1 shows these sets of r and minimum 9 , where 9 is given in degrees. 

Each model is generated by selecting the desirable reduction factor r (r < 1.0), 
the corresponding branch angle 9 given by equation 4.1, the initial length from 
which the reduction will start, the value of N (which in the investigated case is 2), 



31 



and the number of the desired levels for branch generation. Each level contains 2 1 
points corresponding to the end points of the reduced length branches. 



TABLE 4.1 

NUMERICAL VALUES OF r AND MINIMUM ANGLE 9 . 



REDUCTION FACTOR fr) ANGLE 9 

0.53 29.940 

0.55 46.43° 

0.60 95.89° 

0.66 145.35° 

0.71 180.00° 



Figures 4. 1-4.5 show the planar fractal trees that correspond to each set of r 
and 9 of Table 4.1 respectively. As the branch angle 0 increases, the spreading of the 
branches becomes larger. In all trees the physical length of the initial branch is 
chosen as two centimeters. These models were generated by a FORTRAN program 
which has the following input data: 

1. The initial length. This the length of the first dipole whose length will 
be reduced by the constant reduction factor (r). 

2. The value of N, being 2 in all cases. 

3. The initial point that the reduction starts. This point is the end point 
of the initial wire. 
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4. The number of desirable levels for branch generation. 

5. The branch angle 9 . 

The program generates two output files. The first is the input file for the RCS 
program containing all the required data that were mentioned in Chapter 3. The 
second output file contains the coordinates of the start and end points of each 
generated branch. 

For the models characterized by reduction factors 0.53, 0.55, and 0.60, the 
program limits the structure size when the length of a branch becomes smaller than 
0.1 of the wavelength of the incident plane wave. For the models characterized by 
reduction factor 0.66 and 0.71, this limit is taken as 0.15. The reason is that as the 
reduction factor increases the number of branches in a particular branch level 
increases. The size of the tree is truncated so that the number of unknowns is 
manageable. 

The geometry of all planar fractal trees is in the same coordinate system that 
the RCS program uses. All models have been generated in the y^z plane of a 
cartesian coordinate system. 
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Figure 4.1 Fractal Tree for r = 0.53 and 8 = 29.94°. 
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Figure 4.2 Fractal Tree for r = 0.55 and 9 = 46.43°. 
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Figure 4.3 Fractal Tree for r = 0.60 and 0 = 95.89°. 
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Figure 4.4 Fractal Tree for r = 0.66 and 0 = 145.35°. 
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Figure 4.5 Fractal Tree for r = 0.71 and 9 — 180°. 
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V. NUMERICAL RESULTS 



In this chapter, computed results for the radar cross section a of planar 
structures are presented. In order to check the program, the radar cross section 
results of a single planar dipole are compared with the results given in [Ref. 8]. 
Details of this comparison are presented in the following section. 

In all models where numerical results for the radar cross section are presented, 
the following factors have been considered: 

1. E— plane is the plane corresponding to po = 0°, and 6q varying from 0° to 

1800. 

2. H— plane is the plane corresponding to Oq = 90°, and po varying from —90° 
to +90°. 

3. The resulting RCS a is normalized to the square of the wavelength A of the 
incident plane wave (a/ A 2 ). 

4. The number of segments, that each dipole is subdivided, is taken as four 
per wavelength. 

5. The physical dimensions and orientations of the dipoles used to construct a 
planar structure are the same when this structure is investigated at different 
frequencies, and only the segmentation is different depending upon the wavelength A 
of the incident plane wave. 

6. Each PWS mode has length equal to the length of two segments. 

7. The incident electric field is linearly polarized and the parameters a and b 
are taken as 0 and 1 respectively for this investigation. This corresponds to having 
only a horizontal magnetic field. 
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A. SCATTERING FROM A SINGLE DIPOLE 

In Figure 5.1 a centered loaded vertical planar thin dipole of length L and 
width 2w is oriented along the z axis. This dipole is excited by a plane wave of 
frequency 30 GHz traveling in the x— y plane (0o = 90°, ipo = 0°). The same 
situation is described in [Ref. 8, pp. 510—515] for a cylindrical dipole of radius a and 
length L such that L/2a=74.2. In the present case of a planar dipole, the length of 
the planar dipole is selected such that L/2w = 33 [Ref. 9]. Figure 5.2 shows, on a 
semi— log scale, the results computed by the developed program of the monostatic 
normalized radar cross section < 7 / A 2 of this single dipole for values of L/A ranging 
from 0 to 1.4, and for loads Zl = 0 and Zl = 00 . The case of Zl = ® is achieved by 
setting a small gap (0.01 wavelength) at the center of this planar dipole. In Figure 

o 

5.2, the corresponding values of aj A from [Ref. 8, pp. 115] are also shown. 

This comparison leads to an accurate check of the correct calculation of the 
radar cross section by the developed RCS program. The small discrepancy between 
the computed values and those given by [Ref. 8] is due to the error incurred in 
reading values off the curves presented in [Ref. 8, pp. 115]. 
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Figure 5.1 Vertical Centered Loaded Planar Dipole. 
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Figure 5.2 Normalized RCS cr/A^ 



of a Center Loaded Planar Dipole. 
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B. SCATTERING FROM PLANAR FRACTAL TREES 



In Chapter 4, five types of planar fractal trees were described. Each type of 

these trees was lying in the y-z plane. The developed RCS program calculates the 

broadside radar cross section o of each tree for frequencies* 15 GHz— 75 GHz, for the 

monostatic case. The initial dipole that is being reduced by the reduction factor r 

has physical length 2 cm and only the branch lengths are changed for each tree 

depending upon the reduction factor r. As the frequency of the plane wave changes, 

the electrical length of this dipole as well as the dipoles composing the tree will also 

change. For the frequency range of 15 GHz— 75 GHz, the electrical length of the 

initial length will vary from one to four wavelengths. 

o 

Figures 5.3— 5.6 show the variations of ajX , in dB, in the E— plane and the 
H— plane for the case of a fractal tree characterized by reduction factor r = 0.53 and 
branch angle 9 = 29.94°. The frequency of the incident plane wave varied from 15 
GHz to 60 GHz in steps of 15 GHz. The tree is composed of 31 dipoles. This number 
is small as the reduction factor is small and the branch lengths are reduced to very 
small values at low levels. 

2 

The (rfX variations in the E— Plane are in a range of 10 dB approximately. In 
o 

the H-Plane. The ojX is symmetric about the 90° axis and is smoother at lower 

2 

frequencies. As the frequency increases, the maximum value of o\\ at 9q = 90° and 
^o=0°, increases from 2.73 dB at 15 GHz to 19.23 dB at 60 GHz. Figure 5.7 shows 
2 

the variation of the maximum a/X in terms of frequency increments for the same 
structure. This variation is very small between 30 and 45 GHz. 
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Figure 5.3 Normalized RCS in dB at F = 15 GHz of a Fractal Tree 
with r = 0.53 and 6 — 29.94°. 
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Figure 5.4 Normalized RCS at F = 30 GHz of a Fractal Tree 
with r = 0.53 and 8 = 29.940. 
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Figure 5.5 Normalized RCS at F = 45 GIIz of a Fractal Tree 
with r = 0.53 and 0 = 29.94°. 
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Figure 5.6 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.53 and 9 = 29.94°. 



47 



X}- 

n G) 

LO ® 

6 c\j 



CD 




(9P) y / -o 

6 



o 

Figure 5.7 Variation of maximum a/X vs. Frequency 
of a Planar Fractal Tree with r = 0.53 and 6 =29.94°. 
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o 

Figures 5.8—5.12 show the variations of a/ A , in dB, in the E and H planes 

for the fractal tree characterized by a reduction factor r = 0.55 and branch angle 0 

= 46.43°. In this case, the fractal tree is composed of 63 dipoles and spreads more 

than the previous one, and the lengths of the branches are bigger than the previous 

ones (due to a larger reduction factor of 0.55 instead of 0.53). This results in more 
o 

variations for a/ A in both E and H planes. 

This model is investigated at a frequency range of 15—75 GHz. At all 

9 

frequencies, the <r/A“ varies in a range of 10—12 dB approximately. In the H— Plane 

9 

the variation of cr/A is symmetric about the 90° axis. The maximum value of radar 
cross section varies from 2.78 db at 15 GHz to 22.63 dB at 75 GHz. Figure 5.13 

9 

shows the maximum value of aj A“, corresponding to 9q = 90° and ipo = 0° in terms 
of the frequency of the incident plane wave, varying from 15 GHz to 75 GHz. This 
2 

maximum <?/A increases linearly as the frequency increases except the range of 
45-60 GHz where the variation is very small. 
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Figure 5.8 Normalized RCS at F = 15 GHz of a Planar Fractal 
with r = 0.55 and 6 = 46.43°. 
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Figure 5.9 Normalized RCS at F = 30 GHz of a Planar Fractal Tree 
with r = 0.55 and 6 = 46.43°. 
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Figure 5.10 Normalized RCS at F = 45 GHz of a Planar Fractal Tree 
with r = 0.55 and 9 = 46.43°. 



52 



(dB) a / X (dB) 




Figure 5.1 1 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.55 and 0 = 46.43°. 
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Figure 5.12 Normalized RCS at F = 75 GHz of a Planar Fractal Tree 
with r = 0.55 and 9 = 46.43°. 
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Figure 5.13 Variation of Maximum RCS vs. Frequency 
of a Planar Fractal Tree with r = 0.55 and 0 = 46.43°. 
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Figures 5.14—5.18 show the variations of the normalized radar cross section 

2 

a/X of a planar fractal tree characterized by reduction factor r = 0.60 and branch 

angle 8 — 95.89°. This model is composed of 63 dipoles and has more spreading than 

the two models that were previously investigated. The branches have larger physical 

lengths than the ones in other two models. At 15 GHz, the variation of a/ A 2 in the 

E— Plane is very low and similar to the variation in the other two models at the 

same frequency. In the H— Plane this variation is different as the incident plane 

wave strikes more dipoles from —90° to 90°. 

As the frequency increases, the electrical length of the dipoles is also increased, 

and more lobes appear at 30 GHz and higher frequencies in both the E and the H 
o 

planes. The maximum value of a/ X at 8o= 90° and po = 0° varies from -0.23 dB 
at 15 GHz to 21.60 dB at 75 GHz. Figure 5.19 shows the variation of the maximum 

radar cross section in terms of the frequency of the incident plane wave. In a range 

2 

of 15—60 GHz the values of maximum a/X follow a straight line approximately. In 

2 

the range of 60—75 GHz the variation of maximum a/X is small. 

The planar fractal trees characterized by r = 0.66, 8 = 145.35°, and r = 0.71, 0 
= 180°, have large values of branching angles. The large values of reduction factor r 
generate fractal trees whose branches have large physical lengths compared with the 
lengths of the branches of the previously investigated fractal trees. The result is that 
the electrical lengths of these branches are large also at the range of 15—75 GHz, and 
a large number of modes is required to investigate the last two types of fractal trees. 
It was found that the developed RCS program is not able to calculate the radar 
cross section of fractal trees that are characterized by large values of reduction 
factor r and branching angle 8 due to memory restrictions and other numerical 
problems. 
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o 

For this reason the numerical results for oj A are not presented in this thesis for 
these two cases. 

It is seen by comparing Figures 5.7, 5.13, and 5.19 that there is a frequency 
range over which the variation of maximum RCS is rather small. Furthermore, this 
range of frequencies shifts to higher frequencies as the tree structure is spread from r 
= 0.53 to r = 0.60. The maximum RCS of each of these trees varies in a range of 3 
dB approximately at the same frequency. 

Scattering from an actual tree will not exactly follow all these patterns, but it 
is felt that the trends would generally remain the same. This is specially true for the 
frequency behavior. 
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Figure 5.14 Normalized RCS at F = 15 GHz of a Planar Fractal Tree 
with r = 0.60 and 9 — 95.89°. 
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Figure 5.15 Normalized RCS at F = 30 GHz of a Planar Fractal Tree 
with r = 0.60 and 9 — 95.89°. 
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Figure 5.16 Normalized RCS at F = 45 GHz of a Planar Fractal Tree 
with r = 0.60 and 6 = 95.89°. 



60 



(dB) <j / X (dB> 



20 




CfiO 

Figure 5.17 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.60 and 0 = 95.89°. 



61 



30 





-90 -45 0 45 90 

cpo 

Figure 5.18 Normalized RCS at F = 75 GHZ of a Planar Fractal Tree 
with r = 0.60 and 9 = 95.89°. 
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Figure 5.19 Variation of Maximum RCS vs. Frequency 
of a Planar Fractal Tree with r = 0.60 and 6 = 95.89°. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The radar cross section of natural trees is a complicated problem of 
electromagnetic scattering. The real trees are made up of long, intersecting, and 
lossy objects. The geometry of these objects is not easy to set up in a three 
dimensional space. 

In this thesis, the real fractal trees were approximated by planar fractal trees. 
The radar cross section of these trees was calculated by developing a Fortran 
program. The planar fractal trees that were used for this investigation were 
symmetric planar structures composed of perfectly conducting and non— intersecting 
planar dipoles. The geometry of these structures was generated using fractal theory. 

A. CONCLUSIONS 

The radar cross section of the planar fractal trees was calculated using the 
moment method theory. For a given planar structure composed of planar strips and 
illuminated by a plane wave, the program generates the geometry of the PWS 
modes, calculates the impedance between any two of them, and fills the impedance 
matrix [Z]. The voltage on each PWS mode, due to the incident electric field, is also 
calculated, and the voltage column vector [V] is generated. 

The RCS program uses the moment method theory to calculate the current 
distribution on each PWS mode. The radiation equations of electromagnetic theory 
are used to determine the scattered electric field due to the current induced on each 
PWS mode of the structure. The knowledge of both the given incident 
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electric and calculated scattered electric fields leads to the calculation of the radar 
2 

cross section (cr/A ), for the given angles Oo,ipo of the incident electric field. 

The calculated radar cross section of a single centered loaded vertical dipole 
was compared with the values given in [Ref. 8, pp. 115], for a similar dipole. 
Although the investigation was for a planar dipole and in [Ref. 8] the dipole was 
cylindrical, the discrepancy of the results was very small. 

In this thesis, five models of planar fractal trees were investigated. The 
geometry of these models was generated by a given program which was modified to 
generate the input data for the developed RCS program. 

The broadside radar cross section, for the monostatic case, was calculated for 
three of these planar fractal trees, for a frequency range of 15 GHz to 75 GHz in 
steps of 15 GHz. This investigation showed that the radar cross section varied in a 

2 

range of 10 dB in the E— Plane. In the H— Plane, the variation of of A was 

symmetric about the 90° axis and smooth at low frequencies and small branching 

o 

angles. For higher branching angles and reduction factors, more variations of <r/A 
with the frequency were seen in both the E and H planes. 

The maximum RCS at Oo = 90°, <po = 0° showed a small variation over a 
frequency range. This range was different on each investigated tree. The variation of 
the maximum RCS followed a straight line at the other frequencies. In each tree and 
at the same frequency, the maximum RCS varied in the range of 3 dB 
approximately. 

It was found that the developed RCS program was not able to calculate the 
radar cross section of fractal trees that were characterized by large values of 
reduction factor and branching angle due to memory restrictions and other 
numerical problems. 
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The scattering from an actual tree will not exactly follow the patterns that 
were described in Chapter 5, but it is felt that the trends would generally remain 
the same. This is specially true for the frequency behavior. 

B. RECOMMENDATIONS 

One recommendation is to investigate the radar cross section of planar fractal 
trees characterized by values other than those used in this thesis. Especially, a 
limited set of reduction factor and branching angle has to be established for the 
same physical length of the initial dipole. 

Another recommendation is to investigate the radar cross section of planar 
fractal trees characterized by the same values of reduction factor and branching 
angles as those that were used in this thesis but with different physical and 
electrical dimensions of the fractal trees. 

In this thesis, the branches of the fractal trees were assumed to be perfectly 
conducting planar strips. The trees were considered without leaves. The radar cross 
section of fractal trees composed of lossy planar strips and with leaves should be 
investigated. 

Finally, the generation of a three dimensional fractal tree and its radar cross 
section may be investigated. 
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OO OOOOQOOOOOOOOOOOOOOOOOOOOOQOQ 



APPENDIX A 



PROGRAM RCS 



THIS PROGRAM CALCULATES THE RADAR CROSS SECTION OF A 
PLANAR STRUCTURE COMPOSED OF ARBITRARILY ORIENTED 
PLANAR THIN DIPOLES. THE DIPOLES ARE NOT INTERSECTING. 
THE DIPOLES ARE PERFECT CONDUCTORS. 



GEOMETRY IN Y-Z PLANE 

INPUT DATA. 

F = FREQUENCY IN GHZ 
N\Y = NUMBER OF DIPOLES 

A.B = COORDINATES OF INCIDENT ELECTRIC FIELD ON Y AND Z 
AXIS RESPECTIVELY. 

L\Y = HALF LENGTH OF EACH DIPOLE IN CM 

WAY = HALF WIDTH OF EACH DIPOLE IN CM 

SW. TW = COORDINATES OF THE CENTER OF EACH DIPOLE 

ALONG Y AND Z AXIS RESPECTIVELY 

PSIW = ANGLE IN DEGREES BETWEEN THE Z AXIS (REFERENCE 
AND DIRECTION OF CURRENT FLOW ON EACH DIPOLE 
(POSITIVE ANGLES ARE MEASURED CLOCKWISE 
NAY = NUMBER OF SEGMENTS THAT EACH DIPOLE IS 
SUBDIVIDED 

THITAO, PHIO = ANGLES IN DEGREES OF INCIDENT ELECTRIC 
FIELD. 

OUTPUT DATA: 

NORMALIZED RCS (cr/A“) IN dB 

REAL LS(1:230),S(1:230).T(1:230).PSI( 1:230) 

REAL LEN , P G . SU ,T L\ W ID , A .B . M A G ( 1 :2 30 ) 

INTEGER I,NMAX,INDX(230),NT,NP,L,M,N.K,NW 
COMPLEX Z( 1 :230.1 :230), V( 1 :230),CUR(1:230) 

COMMON/ RC1/ NMAX 
COMMON/ RC2/ LSAYI.S,T,PSI 
COMMON/ RC4/ CUR 

OPEN a T NIT=l.FILE=TNPUTl'.FORM=TORMATTED') 

OPEN UNIT=2.FILE='OUTPUTl\FORM=TNFORMATTED') 

OPEN UNIT=3.FILE=TNPUT\FORM=TORMATTED') 

OPEN UNTT=4.FILE='OUTPUT , .FORM='FORMATTED') 
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READ(1,*) F,NW,A,B 
PI = 4.* ATAN(1.) 

KO = PI*F/15. 

RAD = PI/180. 

CALL ZMATR (F,NW,Z) 

NP=230 

NT=NMAX 

CALL CLUDCP (Z,NT,NP,INDX) 

22 READ (3,*,END=11) THITA0,PHI0 
CALL VOLT(F,THITAO,PHIO,A,B,V) 

CALL CLUBSB (Z,NT,NP,INDX,V) 

DO 54 K=1,NT 
CUR(K)=V(K) 

54 CONTINUE 

CALL RADIAT (F,THITAO,PHIO,NMAX,A,B,RCS) 

WRITE(4,*) THITA0,PHI0,RCS 
PRINT*, RCS 
GOTO 22 
11 CLOSE(3) 

STOP 
END 

SUBROUTINE TO COMPUTE THE MUTUAL IMPEDANCE BETWEEN 
THE PWS MODES OF N ARBITRARILY ORIENTED DIPOLES. 

GEOMETRY IN THE Y-Z PLANE 

SUBROUTINE ZMATR (F,NW,Z) 

REAL F, S( 1 :230), T( 1:230), LG, SG,TG,WI(1:230),PSI( 1:230) 

REAL PSIG, PSIW(1:70),LW( 1:70), H,W,DH,DW, WIDTH, LENG 
REAL H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2,LS(1:230) 

REAL WW( 1:70), SW(1:70),TW( 1:70), SS(1:20),TS(1:20) 

REAL PG,WID,LEN,SU,TU,A,B 
INTEGER TEMP,P,L,NS(1:70),NMAX,I,N,NW 
INTEGER FUN,G,J,GB,GR,K 
COMPLEX Z12,Z(1:230, 1:230) 

COMMON/ RC2 / LS,WI,S,T,PSI 
COMMON/ RC1/ NMAX 
COMMON/ JOHN/ LG, NG,SG,TG, PSIG, LENG 
COMMON/ ZPAR/ H,W,DW,DH 
COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
DO 101 I = 1, NW 

READ(1,*) LW(I),WW(I),SW(I),TW(I),PSIW(I),NS(I) 

101 CONTINUE 
CLOSE(l) 

NMAX = 0 
DO 75 1=1, NW 

NMAX = NMAX + (NS (I) - 1) 

75 CONTINUE 
P=1 
GB=0 
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DO 50 I=1,NW 




SG=SW(I) 

TG=TW(I) 

PSIG=PSIW(I) 

CALL GEOM(SS,TS) 
GR=NS(I)+GB— 1 
M = P 

DO 60 J=1,NS(I)— 1 
PSI(M)=PSIG 
LS(M)=LG 
\VI(M)= WIDTH 
S(M)=SS(J) 
T(M)=TS(J) 
PRINT*, S(M),T(M) 
M = M + 1 



GB=GR 

P=GR+1 

50 CONTINUE 

TEMP=NS(1)— 1 
L=1 

G=NMAX 
DO 80 M=1,G 



IF(M.LE.TEMP)GOTO 85 
L=L+1 

TEMP=TEMP+NS(L)— 1 



K=TEMP 
DO 90 N=M,G 

IF(N.LE.K)THEN 

H=LS(M) 

W=WI(M) 

DW=0 

DH=ABS(M— N)*H 

CALL ZSDIP(F,Z12) 

Z(M,N)=Z12 

PRINT*, Z(M,N) 

WRITE(2) Z(M,N) 

ELSE 

H1=LS(M) 

H2=LS(N) 

W1=WI(M) 

W2=WI(N) 

S1=S(M) 

S2=S(N) 

T1=T(M) 
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CONTINUE 
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CONTINUE 
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T2=T(N) 

PSI1=PSI(M) 

PSI2=PSI(N) 

CALL ZPSUR (F,Z12) 

Z(M,N)=Z12 
PRINT*, Z(M,N) 

WRITE(2) Z(M,N) 

ENDIF 

90 CONTINUE 

80 CONTINUE 
DO 24 M=2,G 

DO 26 N=1,M— 1 

Z(M,N)=Z(N,M) 

WRITE(2) Z(M,N) 

26 CONTINUE 

24 CONTINUE 
RETURN 
END 
C 

C SUBROUTINE TO COMPUTE THE MUTUAL/SELF IMPEDANCE 
C BETWEEN TWO DIPOLES. THE DIPOLES ARE ASSUMED TO BE 
C COPLANAR, IDENTICAL AND PARALLEL. THE DIPOLE TO 
C DIPOLE IMPEDANCE IS COMPUTED AS THE SUM OF FOUR 
C MONOPOLE TO MONOPOLE IMPEDANCES. 

C 

C INPUT PARAMETERS: 

C 

C F = Frequency of operation (GHz! 

C H = Half height of the dipole (cm) 

C W = Half width of the dipole (cm) 

C DH = Longitudinal distance between the two dipoles (cm) 

C DW = Transverse distance between the two dipoles (cm) 

C 

SUBROUTINE ZSDIP (F,Z12) 

REAL H, W, DW, DH, F 
COMPLEX Z12, ZT 
COMMON/ ZPAR/ H,W,DW,DH 
ZT = (0.,0.) 

Z12 = (0.,0.) 

CALL ZSMONP (F, H, W, DW, DH, 0, 1, 0, 1, ZT) 

Zi2 = Z12 + ZT 

CALL ZSMONP (F, H, W, DW, DH + H, 0, 1, 1, 0, ZT) 

Z 12 = Z12 + ZT 

CALL ZSMONP (F, H, W, DW, DH - H, 1, 0, 0, 1, ZT) 

Zi2 = Z12 + ZT 

CALL ZSMONP (F, H, W, DW, DH, 1, 0, 1, 0, ZT) 

Z12 = Z12 + ZT 

RETURN 

END 
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C SUBROUTINE TO COMPUTE THE SELF/MUTUAL IMPEDANCE 
C BETWEEN TWO IDENTICAL, COPLANAR MONOPOLES. THE 
C CURRENT IS ASSUMED TO BE CONSTANT IN THE 
C TRANSVERSE DIRECTION. 

C 

C REF: R. JANASWAMY, A SIMPLIFIED EXPRESSION FOR THE 
C SELF/MUTUAL IMPEDANCE BETWEEN TWO COPLANAR 

C AND PARALLEL MONOPOLES, IEEE T-AP, AP-35, 

C No. 10, pp. 1174-1176, October 1987. 

C 

C INPUT PARAMETERS: 

C 

C F = FREQUENCY IN GHz 
C H = LENGTH OF EACH MONOPOLE (cm) 

C W = WIDTH OF EACH MONOPOLE (cm) 

C D = CENTER TO CENTER SPACING BETWEEN THE TWO 
C MONOPOLES IN THE DIRECTION TRANSVERSE TO THE 
C CURRENT FLOW (cm) 

C HH = CENTER TO CENTER SPACING BETWEEN THE TWO 
C MONOPOLES IN THE DIRECTION OF CURRENT FLOW (cm) 

C 

C 111 = TERMINAL CURRENT OF END 1 OF MONOPOLE 1. 

C 121 = TERMINAL CURRENT OF END 2 OF MONOPOLE 1. 

C 112 = TERMINAL CURRENT OF END 1 OF MONOPOLE 2. 

C 122 = TERMINAL CURRENT OF END 2 OF MONOPOLE 2. 

C 

C NOTE: 111, 121, 112, 122 can assume values onlv 0 or 1. 

C ICODE = 0, IF D .LE. 4W 
C ICODE = 1, OTHERWISE 

C With ICODE = 0, the expression provided in the above paper is used. 

C With ICODE = 1, a modified form of the expression provided in the 

C above paper is used. (cf. notes) 

C 

C OUTPUT PARAMETERS: 

C 

C Z12 = COMPLEX IMPEDANCE BETWEEN THE TWO SURFACE 
C MONOPOLES. 

C 

SUBROUTINE ZSMONP (F, H, W, D, HH, 111, 121, 112, 122, Z12) 
REAL F, H, W, D, HH, A, B, PI, K0, V, UB, UBP, UA, UAP 
REAL KW, KH, KD, RC, FR, FI, II, 12, 13, 14, SI, Cl 
REAL UABP, ZR, ZI, X, AA (1), BB (1), SQXV, UAB, TINY 
INTEGER 111, 121, 112, 122, M, N, NX, KI, ICODE 
COMPLEX Z12, AC (-1:1, -1:1), El, E2, Zl, J, CMN, E3, El, FAC 
EXTERNAL FR, FI 

COMMON /PARAM/ N, V, KD, A, B, ICODE 
El (X) = Cl (ABS (X)) — J * SI (X) 

SQXV (X) = SQRT (X * X + V * V) 

TINY = l.E— 6 
PI = 4. * ATAN (1.) 

NX = 1 
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KO = PI * F / 15. 

KW = KO * W 
KH = KO * H 
KD m KO * D 
ICODE = 0 

IF (KD .GT. 4. * KW) ICODE = 1 
KI = 10 
C 

C Note that with this choice of KI accurate results are found in the 
C region 0 < D < 4 W, and for D > 20 W. In between these two regions, 

C accurate results are found with KI = 20. Hence the above choice is 

C good only if this code is used for values of D satisfying the above 
C inequalities. 

C 

FAC = CMPLX (1., 0.) 

A = KD-2. * KW 
B = KD + 2. * KW 
IF (ICODE .EQ. 1) GO TO 2 
AA (1) = A 
BB (1) = B 
GO TO 3 

2 AA (1) = -2. * KW 
BB (1) = 2. * KW 

3 II = (121 * COS (KH) — 111) * 112 

12 = (111 - 121 * COS (KH)) * 122 

13 = 121 - 111 * COS (KH ) * 112 

14 = (121 - 111 * COS (KH)) * 122 
J = CMPLX (0.,1.) 

Z1 = CEXP (J * KH) 

AC (-1,-1) = 12 + II / Zl 
AC (-1, 1) = 12 + II * Z1 
AC 1,-1) = 13-14 * Z1 
AC 1, 1) = 13-14 / Z1 

AC 0, -1) = - (AC (-1, -1) * Zl + AC (1, -1) / Zl) 

AC (0, 1) = - (AC (1, 1) * Zl + AC (-1, 1) / Zl) 

RC = 15. / (2. * SIN (KH) * KW) ** 2 
Z12 = (0., 0.) 

DO 1 M =-l, 1 
DO 1 N = -1, 1, 2 
V = K0 * (HH + M * H) 

IF (ICODE .EQ. 0) GO TO 4 
FAC = CEXP (J * N * V) 

CMN = CMPLX (0.,0.) 

GO TO 5 

4 UA = SQXV (A) + N * V 
UAP = UA - 2. * N * V 
UB = SQXV (B) + N * V 
UBP = UB-2. *N* V 
UAB = SQXV (KD) + N * V 
UABP = UAB - 2. * N * V 
El = El (UB) 
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E2 = (0.,0.) 

IF (ABS (A) .GT. TINY) E2 = El (UA) 

E3 = (0.,0.) 

IF (ABS (KD) .GT. TINY) E3 = El (UAB) 

CMN = 0.5 * (B * B * El + A * A * E2 - 2. * KD * KD * E3 + 
k CEXP (-J * UB) * (1. + J * UBP) + CEXP (-J * UA) * 

k (1. + J * UAP) - 2. * CEXP (- J * UAB) * (1. + J * UABP)) 

5 CALL HABER (NX, AA, BB, FR, KI, ZR) 

CALL HABER NX, AA, BB, FI, KI, ZI) 

Z12 = Z12 + AC (M, N) * (CEXP (J * N * V) * CMN + (ZR + J * ZI) 
k * FAC) 

1 CONTINUE 
Z12 = Z12 * RC 
RETURN 
END 

REAL FUNCTION FR (XI, NX) 

INTEGER N, NX, CODE 

REAL X, V, Tl, KD, A, B, XI (NX), TKW, T2, Cl, SI, TINY 
COMPLEX EL J 

COMMON /PARAM/ N, V, KD, A, B, CODE 
El (X) = Cl (ABS (X)) — J * SI (X) 

TINY = l.E-6 
X = XI (1) 

J = CMPLX (0.,1.) 

IF (CODE .EQ. 1) GO TO 1 
Tl = SQRT (X * X + V * V) 

FR = COS (Tl) 

IF (ABS (V) .GT. TINY) FR = FR * (1. - N * V / Tl) 

IF X .LE. KD) FR = FR * A 
IF (X .GT. KD) FR = -FR * B 
GO TO 2 

1 TKW = 0.5 * (B — A) 

T2 = SQRT ((KD + X) * (KD + X) + V * V) 

FR = REAL (El (T2 + N * V)) 

FR = FR * (TKW - ABS (X)) 

2 END 



REAL FUNCTION FI (XI, NX) 

INTEGER N, NX, CODE 

REAL X, V, Tl, KD, A, B, XI (NX), TKW, T2, Cl, SI, TINY 
COMPLEX El, J 

COMMON /PARAM/ N, V, KD, A, B, CODE 
El (X) = Cl (ABS (X)) - J * SI (X) 

X = XI (1) 

TINY = l.E-6 
J = CMPLX (0., 1.) 

Tl = SQRT (X * X + V * V) 

IF (CODE .EQ. 1) GOTO 1 



73 



oo g noon 



FI = -SIN (Tl) 

IF (ABS (V) .GT. TINY) FI = FI * (1. - N * V / Tl) 
IF X .LE. KD) FI = FI * A 
IF (X .GT. KD) FI = -FI * B 
GO TO 2 

1 TKW = 0.5* (B-A) 

T2 = SQRT ((KD + X) * (KD + X) + V * V) 

FI = AIMAG (El (T2 + N * V)) 

FI = FI * (TKW - ABS (X)) 

2 END 



SUBROUTINE TO COMPUTE THE COORDINATES OF THE PWS 
MODES OF EACH DIPOLE IN THE Y-Z PLANE 



SUBROUTINE GEOM (S,T) 

REAL THI,H1,H2,Y,Z,SG,TG,LG,PSIG 
REAL P,Q,S(1:20),T(1:20),LENG 
INTEGER M,NG,K 

COMMON/ JOHN/ LG,NG,SG,TG,PSIG,LENG 
COSD(X) = COS(X*RAD) 

SIND(X) = SIN(X*RAD) 

PI=4.*ATAN(1.) 

RAD=PI/180. 

THI=90— PSIG 
Hl=COSD(THI) 

H2=SIND(THI) 

Y=SG-(LENG*H1) 

Z=TG-(LENG*H2) 

P=LG*H1 



Q=LG*H2 
S(1)=Y+P 
T(1)=Z+Q 
DO 225 K=2,NG— 1 

S(K)=S(1)+(K— 1)*P 
T(K)=T(1)+(K— 1)*Q 
CONTINUE 
RETURN 
END 



SUBROUTINE ZPSUR (F,Z12) 

REAL HI, Wl, H2, W2, PSI, F, YSTAR, ZSTAR, K0, AI (2), BI (2) 
REAL C (3), D (3), RT, ZT, COT, CSEC, KOS, X, PI, PSI1, PSI2 
REAL SI, Tl, S2, T2, RINTG, IMINTG, RESULT, SIND, COSD,SI,CI 
INTEGER NX, KI 
COMPLEX Z12, J 
EXTERNAL RINTG, IMINTG 

COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, K0 
SIND (X) = SIN (X * PI / 180.) 
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COSD (X) = COS (X x PI / 180.) 

DATA AI, BI / 2 * -1., 2*1./ 

PI = 4. x ATAX (1.) 

PSI = PSI2 — PSI1 

IF (ABS (PSILLE.0.4) PSI= SIGN(1.,PSI)*0.4 
J = CMPLX (0.,1.) 

NX = 2 
KI = 3 

K0 = PI x F / 15. 

C (1) = 1. / SIN (K0 * HI) 

ch) = c (l) 

C (2) = -2. x COS (K0 * HI) * C (1) 

D (1) = 1. / SIN (K0 * H2) 

D (3) = D (1) 

D (2) = -2. x COS (KO x H2) * D (1) 

CSEC = 1. / SIND (PSI) 

KOS = COSD (PSI) 

COT = KOS * CSEC 

YSTAR = (S2-S1) * COSD (PSI1) - (T2-T1) * SIND (PSI1) 
ZSTAR = (S2-S1) x SIND (PSI1) + (T2-T1) * COSD (PSIl) 

RT = YSTAR x CSEC - H2 

ZT = YSTAR x COT - HI - ZSTAR 

CALL HABER (NX, AI, BI, RINTG, KI, RESULT) 

Z12 = RESULT 

CALL HABER (NX, AI. BI, IMINTG, KI, RESULT) 

Z12 = Z12 + J * RESULT 
Z12 = -3.75 x Z12 
RETURN 
END 

REAL FUNCTION SI (XI) 

REAL XI 

REAL AF (4), BF (4), AG (4), BG (4), X, X2, X4, X6, XS, FX, GX, 
t PI, SGX 

DATA AF / 38.027264, 265.187033, 335.67732, 38.102495 / 

DATA BF / 40.021433, 322.624911, 570.23628, 157.105423 / 

DATA AG '/ 42.242855, 302.757865, 352.018498, 21.821899 / 

DATA BG / 48.196927, 482.485984, 1114.978885, 449.690326 / 

SI = 0. 

IF (XI .EQ. 0.) RETURN 
SGN = +1. 

IF (XI .LT. 0.) SGN = -1. 

X = ABS (XI) 

X2 = X x X 

IF (X .GE. 20.) THEN 

FX = 1. / X * (1.-2. / X2) 

GX = 1. / X2 * (1. — 6. / X2) 

GO TO 1 
END IF 
X4 = X2 x X2 
X6 = X4 * X2 
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X8 = X4 * X4 
IF (X .LT. 1.) THEN 

SI = X*(l. - X2 / 18. + X4 / 600. - X6 / 35280. + X8 / 3265920.) 

SI = SI * SGN 
ELSE 

FX = (1. + AF (1) / X2 + AF (2) / X4 + AF (3) / X6 + AF (4) / 
k X8) 

FX = FX / (X * (1. + BF (1) / X2 + BF (2) / X4 + BF (3) / X6 + 
k BF (4) / X8)) 

GX = (1. + AG (1) / X2 + AG (2) / X4 + AG (3) / X6 + AG (4) / 
k X8) 

GX = GX / (X2 * (1. + BG (1) / X2 + BG (2) / X4 + BG (3) / X6 + 
k BG(4) / X8)) 

1 PI = 4. * ATAN (1.) 

X = X - AINT (X / (2. * PI)) * 2. * PI 

SI = SGN * (PI / 2. - FX * COS (X) - GX * SIN (X)) 

END IF 
END 
C 
C 

REAL FUNCTION Cl (X) 

REAL AF (4), BF (4), AG (4), BG (4), X, X2, X4, X6, X8, FX, GX 
REAL PI 

DATA AF / 38.027264, 265.187033, 335.67732, 38.102495 / 

DATA BF / 40.021433, 322.624911, 570.23628, 157.105423 / 

DATA AG / 42.242855, 302.757865, 352.018498, 21.821899 / 

DATA BG / 48.196927, 482.485984, 1114.978885, 449.690326 / 

IF (X .LE. 0.)THEN 

PRINT *, 'Invalid argument for Cl (x)', ' x = ', x 

RETURN 

ELSE 

X2 = X * X 

X4 = X2 * X2 

IF (X .GE. 20.) THEN 

FX = 1. / X * (1. — 2. / X2) 

GX = 1. / X2 * (1.-6. /X2) 

GO TO 1 
END IF 
X6 = X4 * X2 
X8 = X4 * X4 
IF (X .LT. 1.) THEN 

Cl = 0.57721566 + ALOG (X) - X2 * (0.25 - X2 / 96. + X4 / 4320. 
k - X6 / 322560.) 

ELSE 

FX = (1. + AF (1) / X2 + AF (2) / X4 + AF (3) / X6 + AF (4) / 
k X8) 

FX = FX / (X * (1. + BF (1) / X2 + BF (2) / X4 + BF (3) / X6 + 
k BF (4) / X8)) 

GX = (1. + AG (1) / X2 + AG (2) / X4 + AG (3) / X6 + AG (4) / 
k X8) 

GX = GX / (X2 * (1. + BG (1) / X2 + BG (2) / X4 + BG (3) / X6 + 
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t 



PI = 4 




X = X - AINT (X / (2. * PI)) * 2. * PI 
Cl = FX * SIN (X) - GX * COS (X) 
END IF 
END IF 
END 



REAL FUNCTION RINTG (X, NX) 

COMPLEX J, El, TERM1, TERM2, INTGR 

INTEGER NX, M, N, P, Q 

REAL * 4 PSI1, PSI2, Si, Tl, S2, T2 

REAL * 4 C (3), D (3), R, Z, RT, ZT, U, V, CSEC, COT, KOS, Y 
REAL * 4 X (NX), HI, H2, Wl, W2, KO, PZQR, SI, Cl 
,TT, ZM, RN, RMN, TESC, TESD, PI, PZQRP 
COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, KO 
El (Y) = Cl (ABS(Y)) — J * SI (Y) 



INTGR = (0.,0.) 

TESC = ABS (C (2) / C (1)) 

TESD = ABS (D (2) / D (1)) 

TT = ALOG (ABS ((l.+KOS)/(l.-KOS))) 

PI = 4. * ATAN (1.) 

DO 3 M = 1, 3 

IF (M .EQ. 2 .AND. TESC .LE. l.E-6) GO TO 3 
Z = (M— 1) * HI 

ZM = -U * Wl * COT + V * W2 * CSEC + ZT + Z 
TERM2 = (0.,0.) 

DO 2 N = 1, 3 

IF (N .EQ. 2 .AND. TESD .LE. l.E-6) GO TO 2 
R = (N — 1) * H2 

RN = -U * Wl * CSEC + V * W2 * COT + RT + R 
IF (ABS (KO * RN) .LE. 0.5E-1) THEN 

PZQRP = KO * ABS (ZM) - AINT (KO * ABS (ZM) / (2. * PI)) * 



TERM1 = CEXP (-J * PZQRP) * TT 
GO TO 4 
END IF 

IF (ABS (KO * ZM) .LE. 0.5E-1) THEN 

PZQRP = KO * ABS (RN) - AINT (KO * ABS (RN) / (2. * PI)) * 
& 2. * PI 

TERM1 = CEXP (-J * PZQRP) * TT 
GO TO 4 
END IF 

RMN = SQRT (RN * RN + ZM * ZM - 2. * RN * ZM * KOS) 
TERM1 = (0.,0.) 

DO 1 P = -1, 1, 2 
DO 1 Q = -1, 1, 2 
PZQR = P * ZM + Q * RN 



C 




t 



2. * PI 
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PZQRP = PZQR * KO - AINT (PZQR * KO / (2. * PI)) * 2. * PI 
TERM1 = TERM1 + P * Q * CEXP (J * PZQRP) * 
k El (KO * (RMN + PZQR)) 

1 CONTINUE 

4 TERM2 = TERM2 + D (N) * TERM1 

2 CONTINUE 

INTGR = INTGR + C (M) * TERM2 

3 CONTINUE 

RINTG = REAL (INTGR) 

END 

REAL FUNCTION IMINTG (X, NX) 

COMPLEX J, El, TERM1, TERM2, INTGR 

INTEGER NX, M, N, P, Q 

REAL * 4 PSI1, PSI2, SI, Tl, S2, T2 

REAL * 4 C (3), D (3), R, Z, RT, ZT, U, V, CSEC, COT, KOS, Y 
REAL * 4 X (NX), HI, H2, Wl, W2, KO, PZQR, SI, Cl 
k ,TT, ZM, RN, RMN, TESC, TESD, PI, PZQRP 

COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, KO 
El (Y) = Cl (ABS(Y)) — J * SI (Y) 
u = xm 

V = X (2) 

INTGR = (0.,0.) 

TESC = ABS (C (2) / C (1)) 

TESD = ABS (D (2) J D (1)) 

TT = ALOG (ABS ((l.+KOS)/(l.-KOS))) 

PI = 4. * ATAN (1.) 

DO 3 M = 1, 3 

IF (M .EQ. 2 .AND. TESC .LE. l.E-6) GO TO 3 
Z = (M— 1) * HI 

ZM = -U * Wl * COT + V * W2 * CSEC + ZT + Z 
TERM 2 = (0.,0.) 

D0 2N= 1,3 

IF (N .EQ. 2 .AND. TESD .LE. l.E-6) GO TO 2 
R = (N— 1) * H2 

RN = -U * Wl * CSEC + V * W2 * COT + RT + R 
IF (ABS (KO * RN) .LE. 0.5E-1) THEN 

PZQRP = KO * ABS (ZM) - AINT (KO * ABS (ZM) / (2. * PI)) * 
k 2. * PI 

TERM1 = CEXP (-J * PZQRP) * TT 
GOTO 4 
END IF 

IF (ABS (KO * ZM) .LE. .5EM) THEN 

PZQRP = KO * ABS (RN) - AINT (KO * ABS (RN) / (2. * PI)) * 
k 2. * PI 

TERM1 = CEXP (-J * PZQRP) * TT 
GO TO 4 
END IF 

RMN = SQRT (RN * RN + ZM * ZM — 2. * RN * ZM * KOS) 



78 



TERM1 = (0.,0.) 

DO 1 P = -1, 1, 2 
DO 1 Q = -1, 1, 2 
PZQR = P * ZM + Q * RN 

PZQRP = KO x PZQR - AINT (KO * PZQR / (2. * PI)) * 2. * PI 
TERM1 = TERM1 + P * Q * CEXP (J * PZQRP) * 
k El (KO * (RMN + PZQR)) 

1 CONTINUE 

4 TERM2 = TERM2 + D (N) * TERM1 

2 CONTINUE 

INTGR = INTGR + C (M) * TERM2 

3 CONTINUE 

IMINTG = AIMAG (INTGR) 

END 

C 

C Subroutine to compute a sequence of estimates EST1 (K) and 
C EST2 (K), 1 .LE. Kl .LE. K .LE. K2 for the N-dimensional integral 
C B1 BN 

C Int ... Int FUN (xl, x2, ... xN) dxl dx2... dx3 

C A1 AN 

C by Haber's method. 

C Ref: P. J. Davis and P. Rabinowitz, Methods of Numerical 

C Integratation, Academic Press, 1984. 

C 

C For each estimate EST1 (K), two additional quantities ERR1(K) 

C and DEVI (K) are computed. If the values of DEVI (K) do not 
C vary by more than 10% between consecutive values of K, then 
C ERR1 (K) can be taken as a reliable bound on the difference 

C between EST1 and the integral. A similar situation holds for 

C EST2, DEV2, and ERR2 (K). The total number of functional 
C evaluations is 4 * (Kl ** N + (Kl+1) ** N + ... + K2 ** N) and K2 

C should be chosen so as to make this may be halfed by eliminating the 

C computation of the EST2 (K). In other situations, these values 
C are much better than the EST1 (K). A program FUNCTION FUN (X, N) 

C must be supplied by the user with X declared by the statement 

C DIMENSION X (N). FUN must be declared EXTERNAL in the calling 

C program. If N < 1 or N > 10 or Kl < 1 or K2 < Kl, the program 

C terminates with IND = 0. Otherwise IND = 1. 

C 

C Modified by R. Janaswamy so that the output is average of EST1 
C EST2. Also, Kl = K2 = K. 

C 

SUBROUTINE HABER (N, LL, UL, FUN, K, RESULT) 

INTEGER N, IND, KEY, I, K, J 
DOUBLE PRECISION AL (10), BE (10), GA (10), B, G 
REAL FUN, Yl, Y2, Y3, Y4, EST1, EST2, ERR1, ERR2, 
k DEVI, DEV2, RESULT 

REAL * 8 Si, S2, Dl, D2 

REAL LL (N), UL (N), DEX (10), PI (10), P2 (10), P3 (10), P4 (10), 
k Q1 (10), Q2 (10), Q3 (10), Q4 (10), RAN (10), AKN, AK, T, JAC 
REAL AK1, BK 
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EXTERNAL FUN 

DATA AL/. 4142135623730950, .7320508075688773, .236067977499789 
.6457513110645906, .3166247903553998, .6055512754639893, 
.1231056256176605, .3589989435406736, .7958315233127195, 
.3851648071345040/ 

IND = 0 

IF (N .LT. 1 .OR. N .GT. 10) RETURN 
IND = 1 
JAC = 1. 

DO 1 I = 1, N 
BE (I) = AL (I) 

GA (I) = AL (I) 

RAN (I) = UL (I)-LL (I) 

JAC = JAC * RAN (I) 

DEX (I) = 0. 

AK = FLOAT (K) 

KEY = 0 
AK1 = AK-1.1 

51 = 0. 

D1 = 0. 

52 = 0. 

D2 = 0. 

AKN = AK ** N 

T = SQRT (AKN) * AK 

BK = 1. / AK 

KEY = KEY + 1 

IF (KEY .EQ. 1) GO TO 6 

KEY = KEY - 1 



J = 1 

IF (DEX (J) .GT. AK1) GO TO 8 
DEX (J) = DEX (J) + 1. 

GOTO 6 
DEX (J) = 0. 

J = J + 1 



IF (J .LE. N) GO TO 4 
GO TO 3 
DO 7 1= 1,N 
B = BE (I) + AL (I) 

IF (B .GT. 1.) B = B — 1. 

G = GA (I) + B 

IF (G .GT. 1.) G = G -1. 

BE (I) = B + AL (I) 

IF (BE (I) .GT. 1.) BE (I) = BE (I) - 1. 

GA (I) = BE (I) + G 

IF (GA (I) .GT. 1.) GA (I) = GA (I) - 1. 

PI (I) = (DEX (I) + G) * BK 

Ql (I) = LL (I) + RAN (I) * PI (I) 

P2 (I) = (DEX (I) + 1. — G) * BK 
Q2 (I) = LL (I) + RAN (I) * P2 (I) 

P3 (I) = (DEX (I) + GA (I)) * BK 
Q3 (I) = LL (I) + RAN (I) * P3 (I) 
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P4 (I) = (DEX (I) + 1. - GA (I)) * BK 
7 Q4 (I) = LL (I) + RAN (I) * P4 (I) 

Y1 = FUN (Ql, N) 

Y2 = FUN (Q2, N) 

Y3 = FUN (Q3, N 
Y4 = FUN (Q4, N) 

51 = SI + Y1 + Y2 

D1 = D1+ (Y1 - Y2) ** 2 

52 = S2 + Y3 + Y4 

D2 = D2 + (Y1 + Y3 - Y2 - Y4) ** 2 
GO TO 5 

3 EST1 = 0.5 * SI / AKN 

ERR1 = 1.5 * DSQRT (Dl) / AKN 
DEVI = ERR1 * T 
EST2 = 0.25 * (SI + S2) / AKN 
ERR2 = 0.75 * DSQRT (D2) / AKN 
2 DEV2 = ERR2 * T * AK 

RESULT = 0.5 * (EST1 + EST2) * ABS (JAC) 

RETURN 
END 

SUBROUTINE VOLT WHICH CALCULATES THE VOLTAGE MATRIX 
VM 



GEOMETRY IN THE Y-Z PLANE 



SUBROUTINE VOLT(F,THITAO, PHIO, A, B,V) 

REAL TSI,TCO,PS,PC,KY,KZ,A,B,KO,F,KX 
REAL MAR,GTI,PRS,PRC,MOD,KZHTA 

REAL S(1:230),T(1:230),PSI(1:230),LS(1:230),WI(1:230),THITA0,PHI0 

INTEGER NMAX 

COMPLEX V( 1:230), J,STR,VOM 

COMMON/ BRAVO/ KX,KY,KZ 

COMMON/ RC2/ LS,WI,S,T,PSI 

COMMON/ RC1/ NMAX 

COSD(X)=COS(X ,c RAD) 

SIND(X)=SIN(X*RAD) 

PI=4.“ATAN(1.) 

RAD=PI/180. 

K0=PI*F/15. 

J=CMPLX(0.,1.) 

TSI=SIND(THITA0) 

TCO=COSD(THITAO) 

PS=SIND(PHI0) 

PC=COSD(PHI0) 

KY=K0*TSI*PS 
KX=K0*TSrPC 
KZ=K0*TCO 
DO 234 K=1,NMAX 

PRS=SIND(PSI(K)) 
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PRC=COSD(PSI(K)) 

MAR=KY*S(K)+KZ*T(K) 

STR=CEXP(— J*MAR) 

VOM=(STR/SIN(LS(K)*KO))*(A*PRS+B*PRC) 
KZHTA=KY*PRS+KZ*PRC 
MOD=COS(KO*LS(K))— COS(KZHTA*LS(K)) 
GTI=KZHTA**2— K0**2 
IF(ABS(KZHTA— K0)/K0.LT.1.E— 2.0R. 

& .ABS(KZHTA+K0)/K0.LT.1.E-2)THEN 

V(K)=VOM*LS(K)*SIN(KO*LS(K)) 

ELSE 

V(K)=(2*KO/GTI)*VOM*MOD 

ENDIF 

234 CONTINUE 
RETURN 
END 
C 

SUBROUTINE CLUDCP (A, N, NP, INDX) 

INTEGER NP, N, INDX (NP), I, J, K, P, NMAX 
PARAMETER (NMAX = 230) 

COMPLEX A (NP, NP), TEMP, ETA, W (NMAX) 
REAL AAMAX, DUM 
DO 1 K = 1, N— 1 
AAMAX = CABS (A (K, K)) 

P = K 

DO 2 I = K+l, N 
DUM = CABS (A (I, K)) 

IF (DUM .GT. AAMAX) THEN 
AAMAX = DUM 
P = I 
END IF 

2 CONTINUE 
INDX (K) = P 
DO 3 J = 1, N 
TEMP = A (K, J) 

A (K, J) = A (P, J) 

A (P, J) = TEMP 

3 CONTINUE 

DO 4 J = K+l, N 
W (J) = A (K, J) 

4 CONTINUE 

DO 5 I = K+l, N 

ETA = A (I, K) / A (K, K) 

A (I, K) = ETA 

DO 6 J = K+l, N 

A (I, J) = A (I, J) - ETA * W (J) 

6 CONTINUE 

5 CONTINUE 
1 CONTINUE 

RETURN 

END 
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SUBROUTINE CLUBSB (A, N, NP, INDX, X) 

INTEGER N, NP, I, J, INDX (NP), NMAX 
PARAMETER (NMAX = 230) 

COMPLEX A (NP, NP), X (NP), TEMP, Y (NMAX) 

DO PERMUTATIONS ON THE EXCITATION VECTOR USING THE 
INFORMATION ON THE ROW OPERATIONS DONE IN CLUDCP. 



DO 1 K = 1, N— 1 
TEMP = X (K) 

X (K) = X (INDX (K)) 

X (INDX (K)) = TEMP 
CONTINUE 

FORWARD ELIMINATION 



DO 2 I = 1, N 

Y (I) = X (I) 

DO 2 J = 1, 1-1 

Y (I) = Y (I) — A (I, J) * Y (J) 
CONTINUE 



BACK SUBSTITUTION 



DO 3 I = N, 1,-1 
X (I) = Y (I) 



DO 4 J = 1+1, N 

X (I) = X (I) — A (I, J) * X (J) 

CONTINUE 

X (I) = X (I) / A (I, I) 

CONTINUE 

RETURN 

END 



SUBROUTINE RADIAT TO COMPUTE THE RADAR CROSS SECTION 
OF A PLANAR STRUCTURE COMPOSED OF ARBITRARILY 
ORIENTED PLANAR DIPOLES. 

SUBROUTINE RADIAT (F,THITA0,PHI0,NMAX,A,B,RCS) 

REAL RC,THITA0,PHI0,LS(1:230),WI(1:230),S(1:230),T( 1:230), PI 
REAL CR,F,K0,PI,UN,DM,EM,RAD,P,R,A,B,THITA,PHI,RCS,PSI(1:230) 
REAL KX,KY,KZ,KRON,KG 
INTEGER G, NMAX 

COMPLEX J,NTHI0,NPHI0, TEMPI, TEMP2,CUR(1:230) 

COMMON/ BRAVO/ KX,KY,KZ 
COMMON/ RC2/ LS,WI,S,T,PSI 
COMMON/ RC4/ CUR 
COSD(X) = COS(X*RAD) 

SIND(X) = SIN(X*RAD) 
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PI = 4. * ATAN(1.) 
RAD = PI/180. 

KO = PDF/15. 

UN=120.*PI 

THITA =180.- THITAO 

PHI = 180. + PHIO 

J=CMPLX(0.,1.) 

TEMP1=(0.,0.) 

TEMP2=(0.,0.) 

NTHI0=(0.,0.) 

NPHI0=(0.,0.) 

DO 58 G=1,NMAX 
P1=P 



DM=KO*(S(G)*SIND(THITA)*SIND(PHI)+T(G)*COSD(THITA)) 

EM=K0*(SIND(P1)*SIND(THITA)*SIND(PHI)+ 

+C0SD(P1)*SIND(THITA)) 

P=SIND(P1)*C0SD(THITA)*SIND(PHI)-C0SD(P1)*SIND(THITA) 

TEMPl=P*CEXP(rDM)*CUR(G)/SIN(KO*LS(G)) 

TEMP2=CEXP(J*DM)*CUR(G)*SIND(P1)*COSD(PHI)/SIN(KO*LS(G)) 

IF(ABS(EM-K0)/K0.LT.l.E-2.OR.ABS(EM+K0)/K0.LT.l.E-2)THEN 

KG=LS(G)*SIN(KO*LS(G)) 

ELSE 

KG=-2*KO*(COS(EM*LS(G))-COS(KO*LS(G)))/(EM**2-KO**2) 

ENDIF 

NTHI0=TEMP1*KG+NTHI0 

NPHI0=NPHI0+KG*TEMP2 

CONTINUE 

CR=CABS(NTHI0)**2+CABS(NPHI0)**2 

KRON=(A*KY+B*KZ)/KX 

RC=(K0^*2)*(UN**2)*CR/(4.*PI*(ABS(A)**2+ABS(B)**2+ 

+ABS(KRON)**2)) 

IF((RC/((30./F)**2)).LT.l.E— 6)THEN 
RCS = 10.*ALOG10(1.E— 6) 

ELSE 



RCS=10.*ALOG10(RC/((30./F)**2)) 

ENDIF 

RETURN 

END 
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APPENDIX B 



PROGRAM TREE 



C PROGRAM TREE 

C 

C PROGRAM TO GENERATE A PLANAR FRACTAL TREE WITH 
C REDUCTION FACTOR (COND) AND BRANCHING ANGLE THETA 
C WRITTEN BY T.R. NELSON, PhD, UNIVERSITY OF CALIFORNIA 

C SAN DIEGO, LA JOLLA, CA, 90293, AND MODIFIED BY 

C LT. JOHN DEMIRIS TO GENERATE THE INPUT DATA FOR 
C THE RCS PROGRAM. 

C 

C INPUT DATA: 

C 

C L = NUMBER OF BRANCHING LEVELS. 

C N = NUMBER OF BRANCH SEGMENTS. 

C ILEN = INITIAL LENGTH. 

C ISP = INITIAL STARTING POINT. 

C COND = REDUCTION FACTOR. 

C THETA = BRANCHING ANGLE. 

C 

C OUTPUT DATA: 

C 

C IX, IY = COORDINATES OF FIRST AND END POINT OF EACH 
C GENERATED BRANCH. 

C INPUT DATA FOR RCS PROGRAM 

C 

REAL ZZ (1:5,1: 1024), ISP, ILEN, IX, IY,ITH1 

OPEN (UNIT=1, FILE = 'INGEOM',FORM = 'FORMATTED') 

OPEN UNIT=2, FILE = 'OUT', FORM = 'FORMATTED') 

OPEN (UNIT=10,FILE ='INPUTl',FORM='FORMATTED') 
READ(1,*,END=10) L,N, ILEN, ISP, COND, THETA 
10 CLOSE(l) 

ZZ( 1,1) = ISP 
ZZ (4,1) = ILEN 
ZZ (3,1) = 0. 

ZL = ZZ (4,1) 

ZZ (2,1) = 0. 

ZN = N 

PI = 3.14159265 
IX = 0. 

IY = 0.— ILEN 
IX = ZZ(2,1) 

IY = ZZ(1,1) 

DO 100 LOOP =1,L 
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ZLOOP = LOOP 
NPTS = N**ZLOOP 
INC = 0 

NPREV = N**(ZLOOP— 1) 

PREV = NPREV 
ANGLE = (ZN— 1.0)/2.0 
ZINDX = -ANGLE 
DO 200 J =1, NPTS 
ZJ = J 

IF (INC.GE.N) ZINDX=-ANGLE 
IF (INC.GE.N) INC=0 
IW=NPTS— J+l 
IR = PREV— ZJ/ZN+ 1.0 
ZL=ZZ(4,IR) 

T = ZZ(3,IR) 

X = ZZ 1,IR) 

Y = ZZ(2,IR) 

ZL1 = ZL*COND 

IF (ZL1.LT.0.1) GOTO 300 

TH1 = ZLl/33.0 

T1 = T+THETA*ZINDX 

ZZ(3,IW) = T1 

T2 = T1*PI/180.0 

IX = Y 

IY = X 

TEMPX1=IX 

TEMPY1=IY 

IX, IY ARE THE COORDINATES OF THE FIRST POINT OF THE 
BRANCH ALONG THE X,Y AXIS RESPECTIVELY 

WRITE(2,*) IX, IY 

XI = ZLl*COS(T2)+X 
Y1 = ZL1*SIN(T2)+Y 
IX=Y1 

IY=X1 

IX, IY ARE THE COORDINATES OF THE END POINT OF THE 
BRANCH ALONG THE X,Y AXIS RESPECTIVELY 

WRITE(2,*) IX, IY 

TEMPX2=IX 

TEMPY2=IY 

S AND T ARE THE COORDINATES OF THE CENTER OF EACH 
BRANCH 

S = (TEMPX2+TEMPXl)/2. 

T = (TEMPY2+TEMPYl)/2. 

PSI1=ATAN((TEMPX2— TEMPX1)/(TEMPY2— TEMPY1)) 
PSI=PSI1*180./PI 
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ZLHALF=ZLl/2. 

THlHALF=THl/2 

WRITEflO,*) ZLHALF,TH1HALF,S,T,PSI 

ZZ(1,IW)=X1 

ZZ(2,IW)=Y1 

ZZ(4,IW)=ZL1 

ZINDX=ZINDX+1.0 

INC=INC+1 



200 CONTINUE 

100 CONTINUE 
300 STOP 
END 
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